/***************************************************************************\
 * Copyright (C) by Keio University
 * BilateralFilter.cuh created in 01 2013.
 * Mail : fdesorbi@hvrl.ics.keio.ac.jp
 *
 * BilateralFilterTemplate.cuh is part of the HVRL Engine Library.
 *
 * The HVRL Engine Library is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * The HVRL Engine Library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 \***************************************************************************/

#include <hvrl/gpgpu/cuda/Prototypes.hpp>
#include <hvrl/gpgpu/cuda/Tools.cuh>

namespace hvrl {

namespace cuda {

template<typename T>
__global__
void bilateralFilterKernelGray(T* input, T* output, float* kernel,
		const int radius, const float range_coeff, const float range_divider,
		const unsigned int width, const unsigned int height) {

	// Get the corresponding 2D coordinate
	int x = __umul24(blockIdx.x, blockDim.x) + threadIdx.x;
	int y = __umul24(blockIdx.y, blockDim.y) + threadIdx.y;

	if (x < width && x >= 0 && y < height && y >= 0) {
		float sum = 0.0f;
		float t = 0.0f;

		// Get the center gray value
		float center = float(input[y * width + x]);

		for (int i = -radius; i <= radius; i++) {
			if (y + i >= 0 && y + i < height) {

				for (int j = -radius; j <= radius; j++) {
					if (x + j >= 0 && x + j < width) {

						float curPix = float(input[(y + i) * width + x + j]);
						float eulen = (center - curPix) * (center - curPix)
								* range_divider;
						float g = kernel[abs(i) * (radius + 1) + abs(j)]
								* range_coeff * __expf(eulen);
						t += g * curPix;
						sum += g;
					}
				}
			}
		}
		output[y * width + x] = T(t / sum);
	}

}

template<typename T>
void bilateralFilterKernelGray(T *d_input, T *d_output, float* kernel,
		const int& radius, const float& range_coeff, const float& range_divider,
		const unsigned int& d_width, const unsigned int& d_height,
		const unsigned int& block_x, const unsigned int& block_y) {

	dim3 blockSize(block_x, block_y);
	dim3 gridSize((d_width + blockSize.x - 1) / blockSize.x,
			(d_height + blockSize.y - 1) / blockSize.y);
bilateralFilterKernelGray<<< gridSize, blockSize>>>(d_input, d_output, kernel, radius, range_coeff, range_divider, d_width, d_height);
}

template void bilateralFilterKernelGray<unsigned char>(unsigned char *d_input,
		unsigned char *d_output, float* kernel, const int& radius,
		const float& range_coeff, const float& range_divider,
		const unsigned int& d_width, const unsigned int& d_height,
		const unsigned int& block_x, const unsigned int& block_y);
template void bilateralFilterKernelGray<unsigned short>(unsigned short *d_input,
		unsigned short *d_output, float* kernel, const int& radius,
		const float& range_coeff, const float& range_divider,
		const unsigned int& d_width, const unsigned int& d_height,
		const unsigned int& block_x, const unsigned int& block_y);
template void bilateralFilterKernelGray<float>(float *d_input, float *d_output,
		float* kernel, const int& radius, const float& range_coeff,
		const float& range_divider, const unsigned int& d_width,
		const unsigned int& d_height, const unsigned int& block_x,
		const unsigned int& block_y);
template void bilateralFilterKernelGray<int>(int *d_input, int *d_output,
		float* kernel, const int& radius, const float& range_coeff,
		const float& range_divider, const unsigned int& d_width,
		const unsigned int& d_height, const unsigned int& block_x,
		const unsigned int& block_y);

}

}

